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ABSTRACT 

In this paper we propose a fast method for solving wave guide problems. In par- 
ticular, we consider the guide to be inhomogeneous and allow propagation of waves of 
higher order modes. Such problems have been handled successfully for acoustic wave 
propagation problems with single mode and finite length. This paper extends this con- 
cept to electromagnetic wave guides with several modes and infinite length. The 
method is described and results of computations are presented. 
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1. Introduction 

This paper presents a method to accelerate computations involving electromagnetic wave guide 
problems in two dimensions. The wave guide can be finite in length or infinitely long and can have a 
varying index of refraction. The emphasis here is on guides with infinite length because they are 
difficult to handle. In [4] Kriegsmann describes a similar situation and we shall make comparisons with 
his results. However, the method is applicable to more general wave guide problems. For example, 
acoustic wave guide problems (which are treated by Baumeister in [1] and [2]) and optics problems 
involving guided waves such as fiber optics can be treated by this method. Though the discussion here 
pertains to two dimensions, the method could easily be extended to three-dimensional problems with 
cylindrical symmetry. Another advantage of the method is it can be implemented on a personal com- 
puter. 

The tool we use here is a combination of two ideas. The general wave guide problems are posed 
in the frequency domain, i.e., all associated fields have a single angular frequency to and have the time 
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dependency of the form e ia *. The first idea is to solve the corresponding time dependent problem using 
the "limiting Amplitude Principle." The limiting amplitude principle states that the solution of a hyper- 
bolic partial differential equation (which is the case here), with a periodic forcing function, will 
approach a time periodic solution with the same period as t-»°°. This is a known concept; however, a 
rigorous justification in general, including our case, is a difficult task. We refer readers to the paper of 
Morawetz [5] for technical details. Even though the time dependent problem is much harder to treat, it 
has its own advantage. The associated wave equation, which is weakly elliptic, possesses interior 
resonant values for which there are no unique solutions. It is usually difficult to locate such resonant 
frequencies. Considering the problem in the time domain conveniently avoids the situation, and the 
solution passes through these resonant values smoothly yielding the desired solution. This is the advan- 
tage of the limiting amplitude principle. 

The next concept we need here is the wave envelope method. The basis of the method heavily 
relies on the following fact: to compute solutions of time harmonic waves, one usually requires finer 
meshes as the frequency becomes larger and the corresponding computational time increases. However, 
the wave contains an envelope which is the amplitude. Usually, the calculation of the amplitude is the 
goal. For example, in acoustic problems, the root mean square ( r.m.s. ) value of the pressure is 
required rather than the actual time history of the wave. Thus, instead of resolving the actual wave 
oscillations, one may consider a coarse mesh on the envelope as an alternative procedure. A typical 
situation is depicted in figure 1.1. This procedure turns out to be closely connected to the geometric 
optics principle. 

Let us introduce problem under consideration. The problem is posed in the frequency domain as 
follows: Determine u e C 2 (f2), such that 

Au + k 2 n 2 (x,y)u = 0 in Q 
u(x,0) = 0 
u(x, 1 ) = 0 

u(x,y) = e '‘sinny asx^-«>, (k! = V k 2 - n 2 ) , 

u satisfies a radiation condition as x - °° , 


(Pf> -I 
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The physical process here says there is a wave incident in the guide (x of magnitude 

ik x 

e 1 sin rty which impinges through a region R of varying index of refraction and decays at °°. Problem 
(P F ) can be derived from Maxwell’s equations for transverse magnetic field with u being the potential 
of the electric field (see [3]). The homogeneous Dirichlet boundary condition on y = 0 and y = 1 indi- 
cate the guide is insulated. 


Solving (Pp) is the goal of this paper. As we mentioned earlier, first we shall consider the prob- 
lem as a time-dependent problem. The boundary conditions that simulate the behavior at +<*> and -«> 
at finite distances will be considered. Then we will convert into a fast mode with wave envelope tech- 
nique. The discretization will then be considered. Finally, numerical examples will be presented. 


2. The Time-Dependent Problem 

We consider the time-dependent problem corresponding to P F : 

Au - n 2 (x,y)u„ =0 in £2 

u(x,0,t) = 0 
(P T ) ■ u(x,l,t) = 0 

u(x,y,t) = e 1 sm rey, as x ->-» 
u satisfies a radiation condition. 

The goal is to solve for the time-dependent variable u from a state of rest. According to the limiting 
amplitude principle, u(x,y,t) should become u(x,y) e -ikt as t -» °°. A sample calculation to show that 
the time harmonic steady state is achieved is shown in Figure 1.3 (a,b,c,d). In this figure, the time his- 
tory of the real part of the electric field is shown over the indicated span of time. Around t = 20, the 
harmonic steady state is achieved. For illustrative purposes and also to show the stability of the solu- 
tion procedure, the results are indicated up to a time t= 50. The structure of the solution remains 
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unchanged. 

So far we have not specified the boundary condition required as x-»°°. This was omitted 
because we need boundary conditions at finite distances for computational purposes anyway. Both the 
inflow conditions and the radiations condition must be specified at finite distances and should indicate 
the behavior at infinity. We consider the truncated region £2 given as a rectangle £2 = [-L,L] x [0,1], 
where L is typically large (Figure 1.4). 

For L large, the governing equation of (P T ) has constant coefficients (i.e., n(x,y) = 1). Thus, 
near x = -L and x = L the solution using separation of variables may be written as 


- . v iktX - ikt . _ -ik n x - ikt . 

u(x,y,t) = e 1 sin ny + £ R nC sin njty 

(2.1) 

n- 1 


A + ik x - ikt . 

u(x,y,t) = £T n e sin n7ty . 

(2.2) 

n= 1 


k n = V k 2 - n 2 Jt 2 , 

(2.3) 


and R n and T„ are identified as the reflection ans transmission coefficients, respectively. For propagat- 
ing waves, k 2 - n 2 jt 2 > 0. For a given k, the largest number n satisfying the inequality is called the 
’mode’. For example, for k = 3.5 we see that n = 1 and we say that there is one propagating mode. In 
such a case k ip (i=2,3,...) become purely imaginary, and the wave components corresponding to the 
coefficients R ; and Tj (i=2,3,...) will decay. They are called evanescent waves. 

For the relation (2.3) we see that for jt < k < 2jt there will be one mode propagating in the 
guide. Also the solution near x = -L and x = L will be dominated by R t and T^ In fact, we can write 
(2.1) and (2.2) in an approximate sense 

fi(x,y,t) = (e* 1 * + Rje ik,x )e~ ikt sin Tty, (near x= -L), (2.4) 

u(x,y,t) = T^* 1 * ” ‘'“sin jty, (near x = L) . (2.5) 

To prescribe (2.4) and (2.5) as boundary conditions at inflow (x = -L) and outflow (x = L), we need 
to know R, and T, apriori. This information is not available since this is part of the solution. On the 
other hand, we can annihilate Rj and T! using appropriate differential operators. The net result is 
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ICi ^ 

B^u s fi x - — fl t = 2ik i e 1 sin rey, (at x = -L), (2.6) 

ki 

B 2 u ■ u x + — u t = 0, (at x = L). (2.7) 

Thus, (2.6) and (2.7) may be taken as the approximate boundary conditions at the inflow and outflow 
respectively. 


For frequencies in the range jt < k < 4 ji, we see from (2.1), (2.2), and (2.3) that the domi- 
nance of R 2 and T 2 become important. The procedure can be extended to this case by retaining one 
more term in the expansions of (2.1) and (2.2). Thus, in the approximate sense 


u(x,y,t) = (e^sin 7ty+ Rie _iklX sin it y + R 2 e _ik2X sin 2jcy) e _flK (near x = -L) (2.8) 

u(x,y,t) = (Tje^sin try + Tje^sin 2ny) (near x = L) (2.9) 

Analogous to (2.6) and (2.7) we obtain second-order boundary operators as follows: 


* 

_d__ 

k2_a_ 

d ki 

* 

d 

dx 

k dt 

dx k 

dt 

_L + 


d ki_ 

* 

d 

dx 

k dt 

dx + k 

dt 


u = 2k t (kj - k 2 ) e* 1 * 


u= 0, 


ikt 


(x= -L), 


(x = +L). 


( 2 . 10 ) 

( 2 . 11 ) 


For higher values of k, this procedure may be repeated to obtain high-order boundary conditions. 


Details will be found in Kriegsmann [4]. While this procedure remains elegant in extracting boundary 
conditions, their stability with respect to the interior equation is difficult to examine. In fact, the boun- 


dary conditions arising from (2.10) and (2.11) were unstable in our numerical calculations, but reason- 


ably accurate results were obtained with the lower order boundary conditions (2.6) and (2.7). This is in 
contrast to the results that are reported in Kriegsmann [4], 


3. Wave Envelope Technique 

Here we want to make a further change in the dependent variable u of (P T ) so that the new vari- 
able will not oscillate rapidly in the axial direction. For simplicity and to fix ideas, let us consider the 
one-dimensional wave motion governed by 

u tt =u xx . (3.1) 

Let us introduce a new variable v(x,t) such that 
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Substitution of (3.2) in (3.1) yields 


u(x,t) = e 11x v(x,t). 


(3.2) 


v,, = v xx + 2qv x + qV (3.3) 

Let us further specialize to the situation of time harmonic fields with a single angular frequency k, i.e., 


This yields 


v(x,t) = e 41 z(x). 


z" + 2q z' + (T | 2 + k 2 )z = 0, 

whose characteristic roots are 


(3.4) 


(3.5) 


\ = —1) ± ik. (3.6) 

Thus, the general solution is given by 

z(x) = e“ T,x (Ae ik ’‘ + Be _ikx ). (3.7) 

We further specialize (3.2) to the situation that the wave motion is unidirectional moving from left to 

right. Then (3.7) becomes 

z(x) = Ae -(1,_ik)x . (3.8) 

At this point, the choice of q is completely arbitrary. Let q = iri*; then (3.8) further reduces to 

z(x) = Ae _i(11 *~ k)x . (3.9) 

For T|* = k, we see that z(x) = A, which yields the amplitude directly. Basically, we require only two 

mesh points for such a calculation. Of course, this is not the problem we need to solve. However, the 

concept is a valuable tool in understanding the method. 

Let us repeat this process for a problem with variable index of refraction n(x). Here the wave 
motion is described by 

n 2 (x) u tt = u xx . (3.10) 

The process described by (3.2) to (3.5) with q = iq* yields 

z xx + 2iq* z x + (n 2 (x) k 2 - T|* 2 )z = 0. (3.11) 

Let us compare this equation with equation (3.5). In order to capture the amplitude, we set k = q*. 

However, it is no longer feasible in (3.11) to do so, due to the variable nature of n(x). Nevertheless, 


V* 
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for numerical computations, one may slow down the variations of the wave suggested by (3.9) with a 
choice of t|* given by 

•n* = k max n(x). (3.12) 

0£(SL 

In such a situation, the illustration prescribed in figure 1.1 holds. Perhaps there are other ways to 
obtain the best choice of T|* which we do not know at this time. The method is known as the wave 
envelope method. 

Now let us extend this concept to problem (P T ) posed in the truncated domain Q. with boundary 
conditions (2.6) and (2.7). To monitor the solution over a period of length 2it in the axial direction, 
we chose t|* = 27t 'n + . We do not make any choices in the y-direction as the oscillations are rather slow 
in that direction. Thus, our transformation is 


u = e -i2jtT,+x p(x,y,t). (3.13) 

Substituting this relation into (P T ) with boundary conditions (2.6) and (2.7), we obtain the following 
problem: 


Find p e C 2 (£2) such that 


(Pe) 


n 2 (x,y)p u = Ap- Mrtrfp*- (2mi + ) 2 p, 
p(x,0,t) = 0 
p(x,l,t) = 0 


.- + Itj ik,x-ikt+2ltitl + x . 

p x = i2jtiy p + — p t + 2ikie sin jcy 

. k t 

p x = ttrttfp- — p, 


in Q 


(at x = -L) 
(at x = L), 


which is our wave envelope problem. 

Historically, the terminology wave envelope arises because we attempt to resolve the envelope of 


the wave rather that the actual oscillation of the wave. 


4. Numerical Procedure 

The procedure consists of a finite difference method via a variation of the finite area method. We 
refer to Baumeister [1], [2] where we obtained guidance in constructing the procedure. The derivation 
begins from the governing equation for (Pe). which is integrated over a cell (Figure 4.1) 
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J if J _ J_ [ n2 (x.y)P« - Ap + i4;ni + p x + (2icri + ) 2 p] dx dy dt = 0, 

1 2 

where a plus sign (+) in the upper limit of integration indicates evaluation at either the upper or right- 
hand boundary of the integration of the cell and a minus sign (-) in the lower limit indicates evaluation 
at the lower or the left-hand boundary of the cell. In this procedure, the integral involving second 
derivatives in time is approximated by a standard central difference scheme, i.e., 


J. M 2 J_ + J- + -0 dxd y dt = 


Py 1 ~ 2pjj + pfr 


k- 1 


At 


■A, 


( 4 . 2 ) 


Ax 


where A is the area of the cell. For boundary cells 3 and 4, A = -^-Ay, for interior cells 1 and 2, 
A = Ax Ay. 

For an interior point the integration is carried out in a straightforward manner. As in [2] we obtain an 
explicit difference scheme 


D-^ 1 = 
Pi.j 


+ i 


A X. 

Ax 


At 
Ayn 
2 


12 


Ax. 

Ax 


Ax. 

Ax 


+ i2jfn + 
-(2nri + Ay) 2 - 2-2 


- i2jcn + 
21 


Ax. 

Ax 


►Pi+i.j 


f Pi- 1, j + Pi!j+i + Py-t 


AX. 

Ax 


+ 2 


[ALj 

At 


y 


Pu 1 - 


( 4 . 3 ) 


The inflow boundary condition, (see P E ) must be imposed correctly. We shall show this procedure. 
The same procedure holds at the outflow boundary and for other cells. We consider the spatial deriva- 
tive part of (4.1) at the inflow 

Si 
2 


t+ + + r "i 

J m J- J_ [P** + Pry - i4 nn + Px- (2jtti + ) 2 pJdx dy dt. 


( 4 . 4 ) 
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We shall decompose this into four parts: I,, I 2 , I 3 , and I 4 for each integral in the expression 

»+ + + 

r 2 r r ... 


f + 

V 

+ _ 

i£ 

J - 

dx 

. - 


3x 


dy dt 


«+£. + 


=j. 


Pi^ij- p& 


AX 


+ 2i k^* 1 * 4t+ 2 * n| ’'sin Jtyjjdy dt 
_ f t+ 2 JPi+i.j~Py « . v k i dp 1 , 

~K-f ft Ax ~ 2 PiJ_ T"at J Ay 

2iki ik.x - ikt + 2jciT) + x j , . , ,11 

~ — e ^cos(7ty + ) - cos(7ty_) dt 


Pi,, - pfi o + wax) Uj* - 

J Ax k 2 

- jcos(xy.) - cos(xy_)}^- e*'” [e’“* - e M ] 


where 


(4.5) 


y+ = jAy + y Ay, 
y- = jAy - yAy, 
t + = kAt + yAt, 
t_ = kAt - y At, 


and 


Pij = p(iAx,jAy,kAt). 


By a similar calculation we obtain 


t+ Ai 


h= j a. 2 /- J- Pyy dx d y dt 


Pij+i ~ 2p i k j + p k j_i Ax 


Ay 2 

. 4 - Al . . 

2 • + ‘ + 


At, 


I 3 = i4nri + J At J_ p* dx dy dt 


= i47H] + Ay At] 


Pi+U - Py 


2 
+ . + 


I 4 = (2jcti + ) 2 | At J_ J_ p dx dy dt 


(4.6) 


(4.7) 
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= ( 2 jrn + ) 2 P g AyAt. 


Therefore, ( 4 . 4 ) = Ii + I 2 - I3 - I4. 

Repeating the above calculations in cells 1 through 4 , we have 

Pi,j = ^b m p i+ ij + C m Pi_j > j+ + 6 m Pi,j— 1 + fmPi.j j gmPi,j b m 

Where the coefficients am, b m , c m , d m , e m , f m , gm and h m are given in the following table: 


REGION 

a m 

b m 

Cm 

d m 

Cm 

fm 

Sm 


GENERAL 

a l 

bi 

Cl 

1.0 

1.0 

ft 

-1.0 

0 

OBSTACLE 

a 2 

bi 

Cl 

1.0 

1.0 

h 

-1.0 

0 

INLET 

a 3 

b 3 

0 

Ax 

2 Ay 

Ax 

2 Ay 

U 

g 3 

^3 

EXIT 

a 3 

0 

C 4 

Ax 

2 Ay 

Ax 

2 Ay 

h 

g 3 

0 


and where 


a i = 


a 2 = 


a 3 - 


At 

Ayn 

At 

Ayn 

2 At 2 


Ay 

, kl 

ttAx + — At 


to 


b, = 


A L 

Ax 


- 2jirn + -^- 
Ax 


b 3 = - 27 ciri + Ay 

Ax 


Ci = 


Ax. 

Ax 


+ 27ciri + 


Af. 

Ax 


C4 = + 2 jtiTi + Ay 

Ax 


fi = 


U = 


-(2jtifAy) 2 - 2-2 
-( 2 ini + Ay) 2 - 2-2 


' * 
Al 

2 

+ 2 

All 

Ax 

At 

Al 1 

2 

+ 2 

Ayn 

Ax 

h j 

At 


( 4 . 8 ) 


( 4 . 9 ) 
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f,= 


-(2rtlV t ') 2 Ax 2 ^y - — - + (TlAt) 2 AxAy 


and n = tj n(x,y). Finally 


-tj 2 Ax + kj 


g3 = 


At 

co 


■q 2 Ax + kj — 
co 


h 3 = - 2k t e 


i(kj + 2icr| + )x _1_ J 


7t 


cosrt(y+Ay) - costt(y-Ay) 


e ~i<o(t+AO _ e -i<fl(t-At) L 

CO 



2At 

Ay 

t| 2 Ax + ki — 
(0 


The stability requirements are the same as in [1] and translate here into 

Al 


At £ max n(x,y) — 
n 


1 + 


Ax 

V. J 


+ n nA y t| + + (rcAyq + ) 2 
Ax 


With the above scheme, the numerical procedure starts from a state of rest and is driven to time har- 
monic steady state as described in Figure 1.3. 


5. Numerical Results 

In this section we show some results to demonstrate the effectiveness of the method described 
earlier. In numerical computations, the choice of tj and q + is crucial. If q + = 0, there will be no 
envelope transformation, and the solution will be sought with the usual finite difference approach. 
Such choices were mentioned in Section 3. Here, we reinforce this through one example, and the rest 
of the choices follow through the same recipe. In equation (3.12) for the one-dimensional situation, 
we had 

1 ]*= k max n(x) = 2rtri + . (5.1) 

a < x £ L 

The two-dimensional analogue of (5.1) applied to the wave guide problem (P E ) is 

q*= kj max_ n(x,y) = 2itr\ + . (5.2) 

(x.y) E Q 

For instance, for the case of one mode propagating in the guide, we know that the frequency k of the 
incident field must be in the range (re, 2?t). If we consider k = then 


kj = Vk 2 - 7t 2 = “7t . 


( 5 . 3 ) 


- 12 - 


yj 5 

Moreover, if the maximum value of n(x,y) is 1 in £1, then ti + = This procedure is then carried 

out for higher order mode propagation and other frequencies. One must realize this choice is approxi- 
mate; it does not hold exactly even in one dimension unless the index of refraction is a constant. 

Our goal in these calculations is to obtain an efficient procedure that will yield results within a rea- 
sonable time and effort, but with an acceptable range of error. In fact, the code was tested on a 
Masscomp minicomputer, Zenith ZF-151 personal computer, and an IBM AT. Results are obtainable in 
under ten minutes of execution time in all cases although the Masscomp and IBM AT produce results 
two to three times faster than the Zenith. Without the wave envelope technique, about sixty grid 
points in the x-direction were needed to obtain satisfactory results. Comparable results with the wave 
envelope technique required only 13 points in the axial direction. We compare our results with the 
results of Kreigsmann [4] using the same format. 

Our first result is a trivial case where the index of refraction of the media n(x,y) = 1. In this 
case, the wave form should be the same in the axial direction at every point from the inflow boundary 
to the outflow boundary. The value of L was chosen to be 3, so that £1 = [-3,3] x [0,1]. The results 
are presented in Figure 5.1. In this figure the numbers printed at the mesh points are ten times the 
total field. Thus, the high values in the integers correspond to high intensity and the low values to low 
intensity of the electric field. Since our calculations yield p(x,y,t), we used our transformation 

p(x,y,t) = e 12 * 1 ^ v(x,y,t) = e i2,c,1+ ■ ik, u(x,y) 

and calculated the absolute value p to obtain u. Thus, the results are given for u. The frequency of 
oscillations is k = Also, reflection and transmission coefficients at the inflow and the outflow 

boundaries for values of y = .3, .5, and .7 are presented. They give a guideline for monitoring the 
accuracy of the solutions. As we see, the error in the calculation of the reflection coefficients is less 
than 1% and of the transmission coefficients 4.5%. 

The next case we considered is that of a dielectric block imbedded in the guide with index of 
refraction i.e., mu(x,y) = — L- for (x,y), R= [-1,1] x [0,1]. Outside, (in R c ), the index of 

refraction is 1. The steady harmonic state results reported in Figure 1.3 correspond to this situation. 
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The time histories were constructed at the midpoint of the wave guide. Again, Q was chosen to be 
[-3,3] x [0,1], and the frequency k is again -y. The value of q + is 0.559. The results are given in 

Figure 5.2. Also the corresponding reflection and transmission coefficients are given. The exact solu- 
tion can be obtained in this situation as pointed out in [4], The error is almost 3%, and they differ by 
at most 2.5% from the results of [4] with second-order boundary conditions. 


5tc 

Finally, the most interesting case of higher-order modes is presented: the value of k is -y, 

13 

which propagates two modes. A dielectric block with n 2 = — is placed in the region 


R = 4 (x,y) \-l<, \<, 1, 0.2 < y $ 0.8 


r For this case, we do not have exact results or comparis- 


ons, and computations are presented in Figure 5.3. 
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Figure 4.1 



-18- 


NUMERIC REPRESENTATION OF PRESSURE ENVELOPE IN GUIDE 


L 3.00; time=20.00; nxy«= 1.00 

No. of time steps = 421; w= 4.712; etap= 0.559 
Obstacle Boundaries; 

xbl=0 .00; xb2= 0.00; ybl- 0.00; yb2= 0.00 


0000000000000 
3333333333333 
6666666666666 
8888888888899 
101010 9101010 9 9 9101010 
10101010101010 91010101111 
101010 9101010 9 9 9101010 
8888888888899 
6666666666666 
3333333333333 
0000000000000 


REFLECTION, TRANSMISSION coefficients 
point R1 T1 

0.30 0.001 1.055 

0.50 0.001 1.055 

0.70 0.001 1.055 


Figure 5.1. Numeric Representation of Pressure Envelope in Guide 
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L 3. CO; time— 2 C. GO; nxy« 

No. of time steps - 421; 
Obstacle Boundaries; 
xbl— 1.00; xb2« 1.00; ybl— 



C . 58 

w— 4.712; etap- 0.559 


0.00; yb2« 1.00 


0000000000000 
6346310000000 
12 5 712 521000000 
16 7 916 731100000 

19 81119 831100000 

20 91120 842100000 
19 81119 8.31100000 
16 7 916 731100000 
12 5 712 521000000 

6346310000000 

0000000000000 


REFLECTION, TRANSMISSION coefficients 
point R1 T1 
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Figure 5.2 
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Figure 5.3 
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